/* ****************************************************************************
 *
 * Copyright (c) Microsoft Corporation. 
 *
 * This source code is subject to terms and conditions of the Microsoft Public License. A 
 * copy of the license can be found in the License.html file at the root of this distribution. If 
 * you cannot locate the  Microsoft Public License, please send an email to 
 * dlr@microsoft.com. By using this source code in any fashion, you are agreeing to be bound 
 * by the terms of the Microsoft Public License.
 *
 * You must not remove this notice, or any other, from this software.
 *
 * THIS CODE AND INFORMATION ARE PROVIDED "AS IS" WITHOUT WARRANTY OF ANY
 * KIND, EITHER EXPRESSED OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE
 * IMPLIED WARRANTIES OF MERCHANTABILITY AND/OR FITNESS FOR A
 * PARTICULAR PURPOSE.
 * ***************************************************************************/


#include "stdafx.h"
#include <pgl/algotools/RandomGenerator.h>

#ifdef _DEBUG
#undef THIS_FILE
static char THIS_FILE[]=__FILE__;
#define new DEBUG_NEW
#endif



unsigned long	CRandomGenerator::jz;
unsigned long	CRandomGenerator::jsr=123456789;
double			CRandomGenerator::m_dPi=asin( 1.0 ) *2;

long			CRandomGenerator::hz;
unsigned long	CRandomGenerator::iz;
unsigned long	CRandomGenerator::kn[128];
unsigned long	CRandomGenerator::ke[256];
float			CRandomGenerator::wn[128];
float			CRandomGenerator::fn[128];
float			CRandomGenerator::we[256];
float			CRandomGenerator::fe[256];

CRandomGenerator::CRandomGenerator( unsigned long jsreed )
{
	zigset(jsreed);
}


/* nfix() generates variates from the residue when rejection in RNOR occurs. */
float CRandomGenerator::nfix(void) 
{	
	static const float r = 3.442620f; 	/* The starting of the right tail */	
	const double s = 0.290476439455995724186811207;  /* value of one-over-r */

	float x, y;
	for(;;)
	{
		x=hz*wn[iz];
		if(iz==0)
		{	/* iz==0, handle the base strip */
			do
			{	
				x=static_cast<float>(-log(UNI()) * s);   
				/* .2904764 is 1/r */				
				y=static_cast<float>(-log(UNI()));			
			} while(y+y<x*x);
			return (hz>0)? r+x : -r-x;	
		}



		/* iz>0, handle the wedges of other strips */		
		if( fn[iz]+UNI()*(fn[iz-1]-fn[iz]) < exp(-.5*x*x) ) 
			return x;
		/* start all over */		
		hz=SHR3();		
		iz=hz&127;		
		if( static_cast<unsigned int>(abs(hz))<kn[iz] ) 
			return (hz*wn[iz]);	
	}

}



/* efix() generates variates from the residue when rejection in REXP occurs. */

float CRandomGenerator::efix(void)
{	
	float x;
	for(;;)
	{	
		if(iz==0) 
			return static_cast<float>(7.69711-log(UNI()));		/* iz==0 */

		x=jz*we[iz];				
		if( fe[iz]+UNI()*(fe[iz-1]-fe[iz]) < exp(-x) ) 
			return (x);

		/* Start all over again */		
		jz=SHR3();		
		iz=(jz&255);		
		if(jz<ke[iz]) 
			return (jz*we[iz]);	
	}
}



/*--------This procedure sets the seed and creates the tables------*/
void CRandomGenerator::zigset(unsigned long jsrseed) 
{	  
	if (jsrseed == 0) jsrseed = 17235321;

	static const double m1 = 2147483648.0;
	static const double m2 = 4294967296.;

	static const double vn=9.91256303526217e-3;
	static const double ve=3.949659822581572e-3;

	double dn=3.442619855899;
	double de=7.697117470131487;
	double tn(dn);
	double q;      
	double te(de);

	int i;		  
	jsr^=jsrseed;

	/* Set up tables for RNOR */	  
	q=vn/exp(-.5*dn*dn);  
	kn[0]=static_cast<unsigned long>((dn/q)*m1);	  
	kn[1]=0;		  
	wn[0]=static_cast<float>(q/m1);	  
	wn[127]=static_cast<float>(dn/m1);
	fn[0]=1.;	  
	fn[127]=static_cast<float>(exp(-.5*dn*dn));		
	for(i=126;i>=1;i--)      
	{   
		dn=sqrt(-2.*log(vn/dn+exp(-.5*dn*dn)));          
		kn[i+1]=static_cast<unsigned long>((dn/tn)*m1);		  
		tn=dn;          
		fn[i]=static_cast<float>(exp(-.5*dn*dn));          
		wn[i]=static_cast<float>(dn/m1);      
	}

	/* Set up tables for REXP */	 
	q = ve/exp(-de);
	ke[0]=static_cast<unsigned long>((de/q) * m2);	  
	ke[1]=0;
	we[0]=static_cast<float>(q/m2);	  
	we[255]=static_cast<float>(de/m2);
	fe[0]=1.f;	  
	fe[255]=static_cast<float>(exp(-de));		
	for(i=254;i>=1;i--)      
	{   
		de=-log(ve/de+exp(-de));          
		ke[i+1]=static_cast<unsigned long>((de/te)*m2);		  
		te=de;          
		fe[i]=static_cast<float>(exp(-de));          
		we[i]=static_cast<float>(de/m2);      
	}
}
